ניתוח אשכולות לישובים ואזורים סטטיסטיים בישראל, לפי מין וגיל

מבוא

כולנו מניחים את ההנחות המוקדמות לגבי כל מגזר אוכלוסייה בישראל, במיוחד בימים אלו. אך במידה והיינו עוורים למגזר, האם ניתן לקבץ צורות יישוב שונות לפי אופי שונה של יישוב?
התשובה יכולה להינתן די בקלות על ידי ניתוח אשכולות באמצעות שפת התכנות הסטטיסטית R, על ידי ביצוע של הניתוח על נתוני גיל ומגדר של יישובים ואזורים סטטסיטיים.
להלן אציג את הניתוח המוצע על ידי.

שימו לב - הניתוח נהיה טכני קמעה לפעמים. זאת מפני שפוסט זה מציב לעצמו שתי מטרות - גם להציג את יכולתיה של שפת התכנות R לקהל דובר העברית, אך גם לבצע ניתוח נתונים מעניין. במידה ואתם משתעממים ממניפולציה של טבלאות, מציע לדלג על קטעים עם רקע בצבע הזה.

מתודולוגיה, או - תכנית עבודה

מה בתכנית?
1. טעינה וניקוי של הנתונים
2. ניתוח אשכולות
3. יצירת פירמידות גילאים
4. הצגת הנתונים על גבי מפה

שתי הערות טכניות לפני שמתחילים:
אני לא מאמין בהורדת המידע למחשב, אלא שימוש בקבצים זמניים בהורדה ישירה מהמקור. מה שאומר שאם כל הספריות הללו מותקנות לכם על המחשב, אתם יכולים להריץ הכל אצלכם ולראות שאינני מחרטט אתכם.
אני אסיר תודה ליונתן גת שהראה כאן איך ליישר עברית בrmarkdown.

טעינת הספריות הרלוונטיות

לניתוח זה נשתמש בחבילות הבאות:
{tidyverse}, המשמשת לעיבוד ו“שינוי צורה” של נתונים.
{httr}, המשמשת לעבודה עם כתובות אתרים ובקשות מידע מהם.
{readxl}, המשמשת לקריאת קבצי אקסל לתוך R.
{cluster}, המאפשרת ניתוח אשכולות במגוון כלים.
{sf}, המאפשרת קריאה ומניפולציה של קבצים המכילים נתונים מרחביים.
{mapview}, המאפשרת הצגה ויזאולית של נתונים מרחביים.
{leafpop}, המאפשרת הצגה של תמונות כפופאפים למפות.
{rmapshaper}, למניפולציה של ישויות מרחביות.
{randomcoloR}, לבחירת צבעים רנדומליים.

:::
library(tidyverse)
library(httr)
library(readxl)
library(cluster)
library(sf)
library(mapview)
library(leafpop)
library(rmapshaper)
library(htmltools)
library(randomcoloR)

טעינת הנתונים

אנו נסתמך על שני מאגרי נתונים - מאגר נתונים אחד אותו ננתח, ומאגר אחד שיסייע בהצגת הנתונים.
המאגר לניתוח הנתונים נמצא בדף אוכלוסייה ומרכיבי גידול ביישובים ובאזורים סטטיסטיים, 2019-2014, ואנו נשתמש בלוח 5 של שנת 2019.
ניתן לראות כי בקובץ של הלמ“ס מופיעות מספר עמודות רלוונטיות - זכרים, נקבות, ואיחודי א”ס 2019.
טבלאות האקסל המנותחות נוריד את הנתונים, ונביט בטבלת הזכרים:

url1 <- "https://www.cbs.gov.il/he/publications/doclib/2017/population_madaf/population_madaf_2019_7.xlsx"
GET(url1, write_disk(tf <- tempfile(fileext = ".xlsx")))
## Response [https://www.cbs.gov.il/he/publications/doclib/2017/population_madaf/population_madaf_2019_7.xlsx]
##   Date: 2021-03-04 21:54
##   Status: 200
##   Content-Type: application/vnd.openxmlformats-officedocument.spreadsheetml.sheet
##   Size: 1.21 MB
## <ON DISK>  C:\Users\User\AppData\Local\Temp\RtmpSayrHV\file347421d672dc.xlsx
df_male <- read_excel(tf, sheet =4,range = "A12:W2854")
df_female <- read_excel(tf, sheet =5,range = "A12:W2854")
df_unify <- read_excel(tf, sheet =6,range = "A9:D203")
df_male
## # A tibble: 2,842 x 23
##    `סמל צורת יישוב` `סמל יישוב` `שם יישוב` `אזור סטטיסטי` `אוכלוסייה בסוף~ גיל  
##    <chr>                  <dbl> <chr>      <chr>          <chr>            <chr>
##  1  <NA>                     NA <NA>       <NA>           <NA>             4-0  
##  2 "סה\"כ אוכלוסיי~          NA <NA>       <NA>           4537918          4714~
##  3 "סה\"כ אוכלוסיי~          NA <NA>       <NA>           4501377          4650~
##  4 "120"                   3000 ירושלים    סך הכל         467359           60301
##  5 "120"                   3000 ירושלים    111            1444             242  
##  6 "120"                   3000 ירושלים    112            1442             177  
##  7 "120"                   3000 ירושלים    113            1992             444  
##  8 "120"                   3000 ירושלים    114            1952             433  
##  9 "120"                   3000 ירושלים    115            2319             546  
## 10 "120"                   3000 ירושלים    116            3539             824  
## # ... with 2,832 more rows, and 17 more variables: ...7 <chr>, ...8 <chr>,
## #   ...9 <chr>, ...10 <chr>, ...11 <chr>, ...12 <chr>, ...13 <chr>,
## #   ...14 <chr>, ...15 <chr>, ...16 <chr>, ...17 <chr>, ...18 <chr>,
## #   ...19 <chr>, ...20 <chr>, ...21 <chr>, ...22 <chr>, ...23 <chr>

ניתן לראות כי למרות שביקשנו טווח מוגדר (R לא יודע לנחש את מיקום הטבלה כל פעם), חלק משמות העמודות היו בעייתים. לכן יש צורך בסידור הנתונים (ועל הדרך אנחנו גם ננקה אותם).

סידור וניקוי

כפי שניתן לראות, R לא אוהב שמות עמודות בעברית או שמות עמודות כפולות בתאים ממוזגים, ולכן יש צורך בשינוי שמות העמודות.
בנוסף, זרקנו כמה עמודות זבל מהתהליך ועמודות שלא קשורות ליישובים - 3 הראשונות.
כל פעם שהמילה סך הכל הופיעה במספר האזור הסטטיסטי עבור ערים - השמטנו את התצפית על מנת להימנע מכפילות.
ניקינו את הנתונים כך שכל התאים עם ערכי האוכלוסיה יהיו מספריים, ושתאים עם ערך לא ידוע בעקבות ההמרה יהיו שווים 0.
לבסוף, איחדנו בין טבלת הזכרים והנקבות, סינון אזורים בהם סכום הגברים והנשים הוא 0, ומחקנו מלל מהעמודה המצביעה על שם האזור.

# setting workable column name
alt_col_names <- c("type_sym","setl_sym","setl_name","tract","pop2019")
ages <- as.list(as.data.frame(t(df_male[1,])))[[1]]
ages <- ages[!is.na(ages)]
ages <- str_replace(ages,"-","_")
ages <- str_replace(ages,"\\+","_")
ages <- paste0("Age_",ages)
# manual fix
ages[18] <- "Age_up_85"
alt_col_names <- c(alt_col_names,ages)
colnames(df_male) <- alt_col_names
colnames(df_female) <- alt_col_names
# discarding useless rows
df_male <- df_male[-c(1:3),] 
df_female <- df_female[-c(1:3),] 
# filtering sums
df_male <- df_male %>% 
  group_by(setl_sym) %>% 
  mutate(n = n()) %>% 
  ungroup() %>% 
  filter(tract != "סך הכל"|n == 1) %>% 
  select(-n)
df_female <- df_female %>% 
  group_by(setl_sym) %>% 
  mutate(n = n()) %>% 
  ungroup() %>% 
  filter(tract != "סך הכל"|n == 1) %>% 
  select(-n)
# convert relevant columns to numeric and fill NA with zero
df_male <- df_male %>% 
  mutate_at(vars(pop2019:Age_up_85),as.numeric) %>% 
  mutate_at(vars(pop2019:Age_up_85),replace_na,replace = 0)
df_female <- df_female %>% 
  mutate_at(vars(pop2019:Age_up_85),as.numeric) %>% 
  mutate_at(vars(pop2019:Age_up_85),replace_na,replace = 0)
# join males and females, filter zero pop, tract text to NA
df <- df_male %>% 
  left_join(df_female, by = c("type_sym", "setl_sym" ,"setl_name" ,"tract" ),suffix = c(".male",".female")) %>% 
  filter(pop2019.male + pop2019.female != 0) %>% 
  mutate(tract = ifelse(tract == "סך הכל", "1",tract))
df
## # A tibble: 2,753 x 42
##    type_sym setl_sym setl_name tract pop2019.male Age_4_0.male Age_9_5.male
##    <chr>       <dbl> <chr>     <chr>        <dbl>        <dbl>        <dbl>
##  1 120          3000 ירושלים   111           1444          242          205
##  2 120          3000 ירושלים   112           1442          177          133
##  3 120          3000 ירושלים   113           1992          444          352
##  4 120          3000 ירושלים   114           1952          433          351
##  5 120          3000 ירושלים   115           2319          546          472
##  6 120          3000 ירושלים   116           3539          824          648
##  7 120          3000 ירושלים   121           2398          240          199
##  8 120          3000 ירושלים   122           1580          196          137
##  9 120          3000 ירושלים   123           2311          261          188
## 10 120          3000 ירושלים   124           1098           97          101
## # ... with 2,743 more rows, and 35 more variables: Age_14_10.male <dbl>,
## #   Age_19_15.male <dbl>, Age_24_20.male <dbl>, Age_29_25.male <dbl>,
## #   Age_34_30.male <dbl>, Age_39_35.male <dbl>, Age_44_40.male <dbl>,
## #   Age_49_45.male <dbl>, Age_54_50.male <dbl>, Age_59_55.male <dbl>,
## #   Age_64_60.male <dbl>, Age_69_65.male <dbl>, Age_74_70.male <dbl>,
## #   Age_79_75.male <dbl>, Age_84_80.male <dbl>, Age_up_85.male <dbl>,
## #   pop2019.female <dbl>, Age_4_0.female <dbl>, Age_9_5.female <dbl>,
## #   Age_14_10.female <dbl>, Age_19_15.female <dbl>, Age_24_20.female <dbl>,
## #   Age_29_25.female <dbl>, Age_34_30.female <dbl>, Age_39_35.female <dbl>,
## #   Age_44_40.female <dbl>, Age_49_45.female <dbl>, Age_54_50.female <dbl>,
## #   Age_59_55.female <dbl>, Age_64_60.female <dbl>, Age_69_65.female <dbl>,
## #   Age_74_70.female <dbl>, Age_79_75.female <dbl>, Age_84_80.female <dbl>,
## #   Age_up_85.female <dbl>

עיבוד מקדים לנתונים

על מנת לייצר בסיס נתונים להזין אותו לפונקציה שמבצעת את אלגוריתם האשכולות, יש לנרמל את הנתונים עבור כל שכבת גיל, לבטא אותם בצורה יחסית לכלל האוכלוסיה באזור המסוים, ולשמור רק את עמודות הגיל.

df_clus <- df %>% 
  mutate(pop2019 = pop2019.male + pop2019.female) %>% 
  mutate_at(vars(starts_with("Age")),~./pop2019) %>% 
  select(starts_with("Age")) 
df_clus
## # A tibble: 2,753 x 36
##    Age_4_0.male Age_9_5.male Age_14_10.male Age_19_15.male Age_24_20.male
##           <dbl>        <dbl>          <dbl>          <dbl>          <dbl>
##  1       0.0817       0.0692         0.0729         0.0405         0.0419
##  2       0.0615       0.0462         0.0684         0.0531         0.0698
##  3       0.113        0.0899         0.0536         0.0184         0.0232
##  4       0.110        0.0895         0.0553         0.0145         0.0196
##  5       0.118        0.102          0.0515         0.0149         0.0140
##  6       0.117        0.0923         0.0433         0.0135         0.0197
##  7       0.0483       0.0400         0.0436         0.0378         0.0332
##  8       0.0605       0.0423         0.0364         0.0290         0.0318
##  9       0.0540       0.0389         0.0300         0.0259         0.0300
## 10       0.0425       0.0442         0.0469         0.0333         0.0285
## # ... with 2,743 more rows, and 31 more variables: Age_29_25.male <dbl>,
## #   Age_34_30.male <dbl>, Age_39_35.male <dbl>, Age_44_40.male <dbl>,
## #   Age_49_45.male <dbl>, Age_54_50.male <dbl>, Age_59_55.male <dbl>,
## #   Age_64_60.male <dbl>, Age_69_65.male <dbl>, Age_74_70.male <dbl>,
## #   Age_79_75.male <dbl>, Age_84_80.male <dbl>, Age_up_85.male <dbl>,
## #   Age_4_0.female <dbl>, Age_9_5.female <dbl>, Age_14_10.female <dbl>,
## #   Age_19_15.female <dbl>, Age_24_20.female <dbl>, Age_29_25.female <dbl>,
## #   Age_34_30.female <dbl>, Age_39_35.female <dbl>, Age_44_40.female <dbl>,
## #   Age_49_45.female <dbl>, Age_54_50.female <dbl>, Age_59_55.female <dbl>,
## #   Age_64_60.female <dbl>, Age_69_65.female <dbl>, Age_74_70.female <dbl>,
## #   Age_79_75.female <dbl>, Age_84_80.female <dbl>, Age_up_85.female <dbl>

הרצת אלגוריתם לניתוח אשכולות

כעת, כאשר הנתונים מסודרים כיאות, ניתן סוף סוף להריץ את הניתוח הרצוי.
קיימים מספר אלגוריתמים לביצוע הניתוח, אנו נשתמש בניתוח אשכלות היררכי על בסיס חלוקה, היינו - בשלב הראשון כל התצפיות מקובצות לאשכול אחד, ולאט לאט הן מתפרקות לאשכולות שונים, על בסיס המרחק שלהן אחת מהשנייה - ככל שתצפית רחוקה יותר מחברותיה, כך היא תתפרק מהאשכול מוקדם יותר. לבסוף, מקבלים חלוקה של כל התצפיות, לפי סוגים של אזורים המתחברים אחד לשני על סמך מרחק. עוד על שיטת החלוקה - בנספח.

z <- diana(df_clus)
plot(z,  main="plotit", which.plot = 2)

k <- 44
clus_table <- table(cutree(z,k = k)) %>%   data.frame()
clusters <- cutree(z,k = 44)

הניתוח שביצעתי מראה שכאשר מחלקים אותו ל44 חלקים, מקבלים אשכולות “מעניינים”. הגרף הנ"ל נקרא דנדוגרמה, והוא נועד להראות מערכות יחסים היררכיות בין אובייקטים. במקרה שלנו, הוא מראה אינטואיטיבית באיזה מרחק האשכולות נפרדים אחד מהשני.

לא ארחיב כאן בהסבר הטכני, ואשלח את הקוראים לבדוק עוד כאן.

השבת האשכולות אל הטבלה

הפלט של התהליך יוצר לנו וקטור המשייך כל תצפית לאשכול מסוים. נחבר את האשכול חזרה אל הטבלה המקורית, ונבחן את התפלגות האשכולות לפי סוג.

df$clusters <- clusters
clus_table %>% ggplot(aes(Var1,Freq)) + geom_col() + labs(x="מספר אשכול",y="כמות")

df <- df %>% 
  group_by(clusters) %>% 
  mutate(n = n()) %>% 
  ungroup() %>% 
  mutate(clusters = ifelse(n < 20, 99,clusters)) %>% 
  select(-n)

ניתן לראות שאנו מקבלים אשכול אחד שלוקח את רוב הדגימות, אשכול אחד בגודל בינוני, עוד ארבע אשכולות קטנים והשאר זעירים. לכן, בהמשך הניתוח, נתייחס לאשכולות הזעירים כאל אחד, שכן הם מבטא מקרים אקזוטיים, וככל הנראה מדובר ביישובים קטנים מאוד/מוסדות.

יצירת פירמידת גילאים

לאחר שביצענו את הניתוח המסווג את הישובים והאזורים הסטטיסטים לאשכולות, נבצע ויזואליזציה של פירמידות הגילאים באזורים השונים.
הויזואליזציה תתחלק לשניים: ויזואליזציה של האשכולות המייצגים, ויזואליזציה עבור כל אזור.
כאן יתבצע גם תהליך של שיום של האשכולות.

הקוד לשרטוט פירמידת הגילאים מבוסס על התשובה הזאת.

eamp <- df %>% 
  nest(starts_with("Age")) %>% 
  mutate(text = paste(setl_name,"אזור:",tract),
         plo = map2(data,text,~t(.x) %>%
                      as.data.frame() %>%
                      magrittr::extract(1) %>% 
                      rownames_to_column("age_split") %>% 
                      separate(age_split,c("del","to","from_dirty"),sep = "_") %>% 
                      separate(from_dirty,c("from","gender"),sep = "\\.") %>%
                      unite("age_group",from,to, sep = "-",remove = FALSE) %>% 
                      mutate(levels = factor(age_group,levels(fct_reorder(unique(age_group),1:length(unique(age_group)))))) %>% 
                      select(-del) %>%
                      ggplot(aes(x= ifelse(gender == "male", -V1,+V1),y = levels, fill = gender))+
                      geom_col() +
                      lemon::scale_x_symmetric(labels = abs) +
                      labs(x = "Population",y = "age group",title = .y) + 
                      theme(plot.title = element_text(hjust = 0.5))
         ))
facet_plot <- df %>% 
  group_by(clusters) %>% 
  mutate(n = n()) %>% 
  group_by(clusters,n) %>% 
  summarise_at(vars(pop2019.male,pop2019.female, starts_with("Age")), sum) %>% 
  mutate(pop2019 = pop2019.male + pop2019.female) %>% 
  select(-pop2019.male,-pop2019.female) %>% 
  mutate_at(vars(starts_with("Age")),~./pop2019) %>% 
  gather(key,value,-clusters,-pop2019,-n) %>% 
  separate(key,c("del","to","from_dirty"),sep = "_") %>% 
  separate(from_dirty,c("from","gender"),sep = "\\.") %>%
  unite("age_group",from,to, sep = "-") %>% 
  select(-del) %>% 
  mutate(levels = factor(age_group,levels(fct_reorder(unique(age_group),1:length(unique(age_group)))))) %>% 
  filter(clusters != 99) %>% 
  ggplot(aes(x= ifelse(gender == "male", -value,+value),y = levels, fill = gender))+
  geom_col() +
  geom_label(aes(x= -0.05 , y = 14, label = paste("n:",n)),fill = "white") +
  geom_label(aes(x= -0.075 , y = 17, label = paste("tot_pop:",pop2019)),fill = "white",size = 4) +
  lemon::scale_x_symmetric(labels = abs) +
  labs(x = "Fraction of population",y = "age group") +
  facet_wrap(~clusters)
facet_plot  

כפי שניתן לראות, התפלגות הגילאים באשכולות השונים שונה בכל אחד מהאשכולות. ניתן לתת שם לכל אחד מהאשכולות השונים:
אשכול מספר 1 מזכיר מאוד פירדמידות גילאים המאפיינות מדינות עולם מתפתח, ולכן נקרא לאשכול זה - “עולם מתפתח”.
אשכול מספר 2 בולט בכך שהוא מכיל שכבה גדולה של הורים צעירים לילדים רבים. לכן נקרא לו “משפחות צעירות ברוכות ילדים”.
אשכול מספר 3 מזכיר מאוד פירמידות גילאים המאפיינות מדינות עולם מפותח, לכן נקרא לו “מאוזן”.

אשכול מספר 4 מכיל משפחות מעט יותר מבוגרות מאשכול מספר 2 המכילות יחסית פחות ילדים, לכן נקרא לו “משפחות בינוניות”.
אשכול מספר 5 מכיל משפחות מבוגרות יותר מאשכול 4, לכן נקרא לו “משפחות מבוגרות”.
אשכול מספר 6 הוא אחד ההיילייטים של הניתוח. לאחר שגיליתי את קיומו של אשכול זה, הייתי צריך לחשוב מה אני רואה. לאחר ששלחתי מייל ללמ"ס כדי להבהיר מספר דברים, נאמר לי כי:
תשובת הלמ"ס

תשובה זו אוששה את חשדי לגבי מהות אשכול זה. מדובר באשכול המלא בפנימיות לבנים, תוך כדי הימצאות של פירמידת גילאים של עולם מתפתח. לכן לאשכול זה נקרא “ישובים עם ישיבה”.
אשכול מספר 10 בולט בכך שהוא מכיל אוכלוסיה צעירה מרובה ויחסית מעט ילדים, לכן נקרא לו “צעירים”.
אשכול 22 הוא מקרה קיצוני של אשכול 11, לכן נקרא לו “מלא צעירים”. שימו לב שלכל אחד מהאשכולות מצורפים כמות האשכולות ומספר האנשים הכולל המתגורר באשכול.

clusters_names <- c("עולם מתפתח",
                    "משפחות צעירות ברוכות ילדים",
                    "מאוזן",
                    "משפחות בינוניות",
                    "משפחות מבוגרות",
                    "ישובים עם ישיבה",
                    rep(NA,3),
                    "צעירים",
                    rep(NA,11),
                    "ים צעירים")  
facet_plot$data <- facet_plot$data %>% 
  mutate(clusters = clusters_names[clusters])
facet_plot

הורדת מידע גיאוגרפי

כעת, לאחר שסיימנו לנתח את האשכולות השונים, חשבתי שיהיה מעניין לראות אותם על מפה, כדי להבחין בדפוסים מרחביים כאלה ואחרים. בשביל לעשות את זה, נוריד את שכבת האזורים הסטטיסטיים והישובים של הלמ"ס מכאן.

ניקוי מידע גיאוגרפי, וחיבורו ליתר הנתונים

כאשר הורדתי את המידע הגיאוגרפי גיליתי כי הפילוח לפי גילאים היה בתוך השכבה, אם כי לא היה פילוח לפי מגדר. לכן, במידה מסויימת, יתכן שחלק מסידור הנתונים היה מיותר. אם כי, במידה ולא היינו מנתחים את האקסל, לא היינו מוצאים את אשכול מספר 6….

כשניסיתי להנגיש את השכבה בצורה אינטראקטיבית, גיליתי שהיא כבדה מדי, ולכן העברתי אותה תהליך של פישוט.
שוב, מאוד טכני - פישוט שמתבצע על הרבה ישויות שחולקות גבול משותף הוא בעייתי, כי אם הוא מתבצע בנפרד הוא עלול לגרום לחפיפות ולחללים בין הישויות. לכן השתמשתי בחבילה המאפשרת פישוט תוך כדי שמירה על הגבולות בין הישויות. עוד על כך בלינק הקודם, חנונים!
כעת נתחיל בניקוי של השכבה הגאוגרפית. איזורים בהם ערך האוכלוסיה הוא NA, יסוננו החוצה מהשכבה.
בנוסף, נבצע חיבור טבלאי של השכבה לטבלאות הקודמות שיצרנו - נתונים וגרפים של פירמידת הגילאים לכל אזור/ישוב.

shp <- st_cast(shp,"MULTIPOLYGON")
shp <- ms_simplify(shp,keep_shapes = TRUE)
maplot <- shp %>%
  filter(!is.na(Pop_Total),!is.na(STAT11)) %>% 
  mutate(STAT11 = as.character(STAT11)) %>% 
  left_join(df, by =c("SEMEL_YISHUV" = "setl_sym", "STAT11" = "tract")) %>% 
  left_join(eamp %>% select(setl_sym,tract,plo), by =c("SEMEL_YISHUV" = "setl_sym", "STAT11" = "tract")) %>% 
  filter(!is.na(clusters)) %>% 
  mutate(clusters = clusters_names[clusters],
         clusters = ifelse(is.na(clusters), "אקזוטי",clusters),
         plo = map2(plo,clusters,~.x + labs(subtitle = .y) + theme(plot.subtitle = element_text(hjust = 0.5))))

הצגת הנתונים על גבי מפה

ועכשיו לרגע שכולנו חיכינו לו - המפה! המפה מראה את גבולות האיזורים ומסמנת לפי צבע את האשכול אליו הם שייכים. בנוסף, לחיצה על איזור מאפשרת לראות את פירמידת הגילאים שלו. הערה קטנה לגבי התנחלויות: למעט הערים בהתנחלויות, הלמ"ס לא שומר על שכבות של הגבולות של כל אחד מהישובים, אלא מסמן אותם כנקודה שעליה עושים חיץ קטן. אם חשקה נפשכם לראות נתונים של התנחלות כזו או אחרת, פשוט תתמרכזו אליה.
מזמין אתכם לגלות דברים מעניינים. למשל - האם אשכולות מסוימים נוטים להיות קרובים אחד לשני? האם אתם יכולים להזדהות עם הניתוח לגבי אזורים מסוימים בעיר שלכם? האם המפה מחדשת לכם דברים שלא ידעתם?
בכל מקרה, חובה עליכם לעשות זום אין על אזור תל אביב. מאוד מעניין מה שקורה שם.

mapviewOptions(vector.palette = randomColor(9))
map <- mapview(maplot,zcol = "clusters",popup = popupGraph(maplot$plo))
map@map

נספח: אלגוריתם לקביעת חלוקה לקבוצות

על מנת לחלק את האשכולות לקבוצות הגיוניות, חיפשתי את הנקודה בה אשכול גדול מתפרק לאשכולות קאנים. נקודות אלו יכול להיות מזוהות על ידי ירידה קיצונית בשונות בין שלבים. הגרף שאני מציג כאן מראה באלו שלבים היה הגיוני לפרק לאשכול שונים. ניתן לראות כי היה אפשר לפרק את הגרף לאשכולות שונים כבר כאשר היו 16 אשכולות, אך אם היינו עושים זאת לא היינו מגלים הרבה דברים מעניינים.

lapply(2:100,function(x){
  table(cutree(z,k = x)) %>%   data.frame(n=x)
}) %>% 
  bind_rows() %>% 
  group_by(n) %>% 
  summarise(ksd = sd(Freq),kmean = mean(Freq),kcv = ksd/kmean) %>% 
  ungroup() %>% 
  mutate(diff = kcv - lag(kcv), rel = ksd / lag(ksd)) %>% 
  gather(metric, value,-n) %>% 
  mutate(sep_facet = ifelse(metric == "ksd" | metric == "kmean",metric,"unify")) %>% 
  ggplot(aes(x=n,y=value,color = metric)) + 
  geom_line() + 
  facet_grid(sep_facet~.,scales = "free_y")
## `summarise()` ungrouping output (override with `.groups` argument)